R Code for Hoefer et al. 2024 - Sensors vs Surveyors: Comparing passive acoustic monitoring, camera trapping and observer-based monitoring for terrestrial mammals.

Author

Sebastian Hoefer

Published

August 13, 2025

Load packages

Load packages needed for data manipulation and analyses.

library(tidyverse)    # Data wrangling and visualization
library(colorspace)   # Color manipulation
library(cmdstanr)     # Interface to 'CmdStan'
library(brms)         # Bayesian regression models
library(rstan)        # Interface to Stan
library(ggeffects)    # Marginal effects plots
library(DHARMa)       # Residual diagnostics
library(emmeans)      # Estimated marginal means
library(tidybayes)    # Tidy data for Bayesian models
library(vegan)        # Community ecology analysis
library(EcolUtils)    # Ecology utilities
library(patchwork)    # Combine ggplots
library(HDInterval)   # HPD intervals
library(report)       # Automated reporting
library(ggsignif)     # Significance levels in ggplot
library(ggimage)      # Images in ggplot
library(cowplot)      # Plot composition in ggplot
library(paletteer)    # For colour palettes
library(pilot)        # For colour palettes

Functions

Custom ggplot theme for visualisation.

my.theme <- function(){
  theme_classic() +
    theme(text = element_text(family = "Avenir Next"),
          axis.title.y = element_text(margin = margin(t = 0,r = 20,b = 0,l = 0)),
          axis.title.x = element_text(margin = margin(t = 20,r = 0,b = 0,l = 0)), 
          plot.margin = unit(c(5, 10, 5, 10), units = "mm"),
          strip.background = element_rect(fill = "#CCCCFF"),
          strip.text.x = element_text(size = 20),
          axis.title = element_text(size = 20),
          axis.text = element_text(size = 18),
          legend.text = element_text(size = 15),
          legend.title = element_text(size = 15))
}

Load data

All mammal data

mammal_data <- read.csv("mammal_data.csv")
mammal_data_vocal <- mammal_data %>% 
  filter(vocal == "yes",
         !common.name %in% c("Yellow-bellied Sheathtail-bat")) # Removed because mostly out of human audible range

Function to summarise the total number of species per method for each plot across all sites. Adding missing zero values for methods that didn’t capture any species at a given plot.

summarise_mammal_data <- function(df) {
  
  site_plots <- c("Tarcutta.DryA", "Tarcutta.DryB", "Tarcutta.WetA", "Tarcutta.WetB", 
                  "Duval.DryA", "Duval.DryB", "Duval.WetA", "Duval.WetB", 
                  "Mourachan.DryA", "Mourachan.DryB", "Mourachan.WetA", "Mourachan.WetB",
                  "Wambiana.DryA", "Wambiana.DryB", "Wambiana.WetA", "Wambiana.WetB",
                  "Undara.DryA", "Undara.DryB", "Undara.WetA", "Undara.WetB",
                  "Rinyirru.DryA", "Rinyirru.DryB", "Rinyirru.WetA", "Rinyirru.WetB")

  sites <- c("Tarcutta", "Duval", "Mourachan", "Wambiana", "Undara", "Rinyirru")
  
  methods <- c("PAM.all", "OBM", "PAM.survey", "camera")
  
  df <- df %>%
    unite(site.plot, c(site, plot), sep = ".", remove = FALSE) %>%
    select(-plot) %>%
    group_by(site, site.plot, assessment.method) %>%
    summarise(richness = length(unique(common.name))) %>% # Needs to be common.name because of Dingo and Dog having the same scientific.name
    ungroup() %>%
    select(site, site.plot, assessment.method, richness) %>%
    full_join(crossing(site.plot = unique(.$site.plot),
                       assessment.method = unique(.$assessment.method)),
              by = c("site.plot", "assessment.method")) %>%
    replace_na(list(richness = 0)) %>%
    mutate(site.plot = factor(site.plot, levels = site_plots),
           site = str_extract(site.plot, "^[^.]*"),
           site = factor(site, levels = sites),
           assessment.method = factor(assessment.method, levels = methods),
           richness = as.numeric(richness))

  return(df)
}

Run function.

mammal_data_summary <- summarise_mammal_data(mammal_data)
mammal_data_vocal_summary <- summarise_mammal_data(mammal_data_vocal)

PAM data

pam_data <- read_csv("mammals_PAM_daily_all.csv", show_col_types = FALSE) %>% 
  filter(MANUAL.ID != "UNSURE") %>% 
  mutate(MANUAL.ID = ifelse(MANUAL.ID == TRUE, 1, 0)) %>% 
  group_by(template) %>%
  filter(any(MANUAL.ID == 1)) %>% #Filter for templates with at least one true positive
  ungroup()

BirdNET Embeddings Performance Check

Custom function

Creating a smaller data frame specifically for plotting a smooth line (100 points per template) in plots.

# Function to return predictions over extended range
run_model <- function(df) {
  model_logistic <- glm(MANUAL.ID ~ eucledian.distance, 
                       family = binomial(link = "logit"), 
                       data = df)
  
  # Create extended x range
  new_x <- data.frame(
    eucledian.distance = seq(min(df$eucledian.distance), 
                           max(df$eucledian.distance), 
                           length.out = 100)
  )
  
  # Get predictions and SE
  preds <- predict(model_logistic, newdata = new_x, type = "link", se.fit = TRUE)
  
  # Add predictions to new_x
  new_x$fit <- plogis(preds$fit)
  new_x$lower <- plogis(preds$fit - 1.96 * preds$se.fit)
  new_x$upper <- plogis(preds$fit + 1.96 * preds$se.fit)
  
  # Add template column to predictions
  new_x$template <- unique(df$template)
  
  return(new_x)  # Return only the predictions dataframe
}


# Get the predictions
model_predictions <- pam_data %>%
  group_by(template) %>%
  group_modify(~run_model(.x)) %>%
  ungroup()

Recall calculations Embeddings

Sample from actual values (does not assume normality).

sample.vec <- function(x, size = 1) {
  x[sample.int(length(x), size = size, replace = TRUE)]
}

# Function that uses a simple for-loop to compare random draws from positives vs negatives
calculate_probability_loop <- function(distancesPositives, distancesNegatives, nComparisons = 10000) {
  
  set.seed(123)
  # Set up a vector to record comparisons
  comparisons <- c()
  for (i in seq_len(nComparisons)) {
    # Draw one random value from positives and one from negatives
    comparisons <- c(
      comparisons,
      sample.vec(distancesPositives, 1) < sample.vec(distancesNegatives, 1)
    )
  }
  # Return the fraction of times pos < neg
  sum(comparisons) / length(comparisons)
}

# Now apply this by template
template_probabilities <- pam_data %>%
  group_by(template) %>%
  summarise(
    probability = calculate_probability_loop(
      eucledian.distance[MANUAL.ID == 1],
      eucledian.distance[MANUAL.ID == 0],
      10000
    )
  )

Visualisation

(combined_plot <- ggplot(pam_data, aes(x = eucledian.distance, y = MANUAL.ID)) +  
  # Add CI ribbon
  geom_ribbon(data = model_predictions,
              aes(y = fit, ymin = lower, ymax = upper),
              fill = "grey70",
              alpha = 0.3) +
  # Add CI boundary lines
  geom_line(data = model_predictions,
            aes(y = lower),
            linetype = "dashed",
            color = "black",
            linewidth = 0.8) +
  geom_line(data = model_predictions,
            aes(y = upper),
            linetype = "dashed",
            color = "black",
            linewidth = 0.8) +
  # Add mean prediction line
  geom_line(data = model_predictions,
            aes(y = fit),
            color = "black",
            linewidth = 2) +
  geom_rug(data = subset(pam_data, MANUAL.ID == 0),
           aes(x = eucledian.distance),
           color = "#5093CB", alpha = 0.3,
           length = unit(0.05, "npc"),
           linewidth = 0.8, 
           sides = "b",
           inherit.aes = FALSE) +
  geom_rug(data = subset(pam_data, MANUAL.ID == 1),
           aes(x = eucledian.distance),
           color = "#F78303", alpha = 0.3,
           length = unit(0.05, "npc"),
           linewidth = 0.8, 
           sides = "t",
           inherit.aes = FALSE) +
  geom_text(data = template_probabilities,
            aes(x = Inf, y = Inf, label = sprintf("%.2f", probability)),
            hjust = 1.3, vjust = 2, inherit.aes = FALSE) +
  facet_wrap(~ template, scales = "free_x", ncol = 10) +
  labs(x = "Euclidean Distance", 
       y = "Probability of correct prediction") +
  scale_y_continuous(limits = c(-0.1, 1.1), breaks = seq(0, 1, 0.25)) +
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.text = element_text(face = "bold"),
        axis.title = element_text(size = 23, face = "bold"),
        axis.title.x = element_text(margin = margin(t = 25, r = 0, b = 0, l = 0)),
        axis.title.y = element_text(margin = margin(t = 0, r = 25, b = 0, l = 0)),
        plot.margin = unit(c(5, 10, 5, 10), units = "mm")))

Mammal Traits Effect

Imported from COMBINE (https://doi.org/10.1002/ecy.3344) and added five species not in the database using Field Companion to the Mammals of Australia and HomeRangeData (https://doi.org/10.1111/geb.13625). HomeRangeData was also used to fill in missing home range information. Detectability was calculated by divided true positive detections by overall detections and added.

trait_detect <- read_csv("trait_detect.csv") %>% 
  mutate(fossoriality = as.factor(fossoriality),
         activity_cycle = as.factor(activity_cycle),
         terrestrial_volant = as.factor(terrestrial_volant))
Rows: 19 Columns: 8
── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): scientific.name
dbl (7): adult_mass_g, dispersal_km, fossoriality, home_range_km2, activity_cycle, terrestrial_volant, detectability

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
custom_labels <- function(x) {
  ifelse(x >= 1, 
         scales::comma(x, accuracy = 1),                    # Whole numbers for ≥ 1
         ifelse(x >= 0.1,
                scales::comma(x, accuracy = 0.1),           # 1 decimal for 0.1 to < 1  
                scales::comma(x, accuracy = 0.01)))         # 2 decimals for < 0.1
}

p1 <- trait_detect %>%
ggplot(aes(adult_mass_g/1000, detectability)) +
  geom_smooth(method = "lm", col = "#F78303", linewidth = 1.5) +
  geom_point(size = 2.5, alpha = 0.8) +
  my.theme() +
  scale_y_continuous(name = "Detection probability",
                     breaks = seq(0, 0.10, by = 0.05)) +
  coord_cartesian(ylim = c(-0.02, 0.12)) +
  scale_x_log10(name = "Adult mass (kg)", 
                labels = custom_labels,
                expand = expansion(mult = c(-0.14, 0.05)))
p2 <- ggplot(trait_detect, aes(dispersal_km, detectability)) +
  geom_smooth(method = "lm", col = "#F78303", linewidth = 1.5) +
  geom_point(size = 2.5, alpha = 0.8) +
  my.theme() +
  scale_y_continuous(name = "Detection probability",
                     breaks = seq(0, 0.10, by = 0.05)) +
  coord_cartesian(ylim = c(-0.02, 0.12)) +
  scale_x_continuous(name = "Dispersal (km)",
                     breaks = seq(0,16,by=2),
                     expand = expansion(mult = c(0.02, 0.05)))
p3 <- trait_detect %>% 
ggplot(aes(home_range_km2, detectability)) +
  geom_smooth(method = "lm", col = "#F78303", linewidth = 1.5) +
  geom_point(size = 2.5, alpha = 0.8) +
  my.theme() +
  scale_y_continuous(name = "Detection probability",
                     breaks = seq(0, 0.10, by = 0.05)) +
  coord_cartesian(ylim = c(-0.02, 0.12)) +
  scale_x_log10(name = "Home range (km)", 
                labels = custom_labels,
                expand = expansion(mult = c(-0.12, 0.05)))
pal <- c("#204366","#CB8D79")
pal_dark <- after_scale(darken(pal, .1, space = "HLS"))

p4 <- trait_detect %>% 
  drop_na(fossoriality) %>% 
ggplot(aes(fossoriality, detectability, fill = fossoriality)) +
  geom_boxplot(aes(color = fossoriality,
                   fill = after_scale(desaturate(lighten(color, .8), .4))),
               position = position_dodge2(width = 0.4),
               width = .6, outlier.shape = NA) +
  geom_point(aes(fill = fossoriality), color = "transparent",
             position = position_dodge2(width = 0.4),
             shape = 21, stroke = .4, size = 3.5, alpha = .4) +
  my.theme() +
  scale_fill_manual(values = pal, guide = "none") +
  scale_colour_manual(values = pal_dark,guide = "none") +
  scale_y_continuous(name = "Detection probability",
                     breaks = seq(0, 0.10, by = 0.05)) +
  scale_x_discrete(name = "Fossoriality",
                   labels = c("Fossorial", "Above ground dwelling"))
pal2 <- c("#204366","#CB8D79", "#324742")
pal2_dark <- after_scale(darken(pal2, .1, space = "HLS"))

p5 <- trait_detect %>% 
  drop_na(activity_cycle) %>% 
ggplot(aes(activity_cycle, detectability, fill = activity_cycle)) +
  geom_boxplot(aes(color = activity_cycle,
                   fill = after_scale(desaturate(lighten(color, .8), .4))),
               position = position_dodge2(width = 0.4),
               width = .6, outlier.shape = NA) +
  geom_point(aes(fill = activity_cycle), color = "transparent",
             position = position_dodge2(width = 0.4),
             shape = 21, stroke = .4, size = 3.5, alpha = .4) +
  my.theme() +
  scale_fill_manual(values = pal2, guide = "none") +
  scale_colour_manual(values = pal2_dark,guide = "none") +
  scale_y_continuous(name = "Detection probability",
                     breaks = seq(0, 0.10, by = 0.05)) +
  scale_x_discrete(name = "Activity cycle",
                   labels = c("Nocturnal", "Cathemeral", "Diurnal"))
pal3 <- c("#204366","#CB8D79")
pal3_dark <- after_scale(darken(pal3, .1, space = "HLS"))

p6 <- trait_detect %>% 
  drop_na(terrestrial_volant) %>% 
ggplot(aes(terrestrial_volant, detectability, fill = terrestrial_volant)) +
  geom_boxplot(aes(color = terrestrial_volant,
                   fill = after_scale(desaturate(lighten(color, .8), .4))),
               position = position_dodge2(width = 0.4),
               width = .6, outlier.shape = NA) +
  geom_point(aes(fill = terrestrial_volant), color = "transparent",
             position = position_dodge2(width = 0.4),
             shape = 21, stroke = .4, size = 3.5, alpha = .4) +
  my.theme() +
  scale_fill_manual(values = pal3, guide = "none") +
  scale_colour_manual(values = pal3_dark,guide = "none") +
  scale_y_continuous(name = "Detection probability",
                     breaks = seq(0, 0.10, by = 0.05)) +
  scale_x_discrete(name = "Volant",
                   labels = c("Non-flying", "Flying"))
(traits_plot <- p1 + p2 + p3 + p4 + p5 + p6 +
    plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 20, face = "bold")))

Investigation

trait_detect_clean <- trait_detect %>% 
  drop_na(detectability)

model1 <- glm(detectability ~ adult_mass_g, data = trait_detect_clean)
summary(model1)

Call:
glm(formula = detectability ~ adult_mass_g, data = trait_detect_clean)

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.876e-02  8.854e-03   4.378  0.00041 ***
adult_mass_g -1.834e-08  3.461e-08  -0.530  0.60299    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.001249275)

    Null deviance: 0.021589  on 18  degrees of freedom
Residual deviance: 0.021238  on 17  degrees of freedom
AIC: -69.212

Number of Fisher Scoring iterations: 2
model2 <- glm(detectability ~ dispersal_km, data = trait_detect_clean)
summary(model2)

Call:
glm(formula = detectability ~ dispersal_km, data = trait_detect_clean)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.039930   0.017131   2.331   0.0481 *
dispersal_km -0.001082   0.002359  -0.459   0.6587  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.001204296)

    Null deviance: 0.0098876  on 9  degrees of freedom
Residual deviance: 0.0096344  on 8  degrees of freedom
  (9 observations deleted due to missingness)
AIC: -35.071

Number of Fisher Scoring iterations: 2
model3 <- glm(detectability ~ fossoriality, data = trait_detect_clean)
summary(model3)

Call:
glm(formula = detectability ~ fossoriality, data = trait_detect_clean)

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)    0.039218   0.010709   3.662  0.00193 **
fossoriality2 -0.005552   0.016504  -0.336  0.74067   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.001261516)

    Null deviance: 0.021589  on 18  degrees of freedom
Residual deviance: 0.021446  on 17  degrees of freedom
AIC: -69.027

Number of Fisher Scoring iterations: 2
model4 <- glm(detectability ~ home_range_km2, data = trait_detect_clean)
summary(model4)

Call:
glm(formula = detectability ~ home_range_km2, data = trait_detect_clean)

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)   
(Intercept)     3.722e-02  9.998e-03   3.723  0.00256 **
home_range_km2 -1.391e-06  5.429e-06  -0.256  0.80183   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.001394961)

    Null deviance: 0.018226  on 14  degrees of freedom
Residual deviance: 0.018134  on 13  degrees of freedom
  (4 observations deleted due to missingness)
AIC: -52.202

Number of Fisher Scoring iterations: 2
model5 <- glm(detectability ~ activity_cycle, data = trait_detect_clean)
summary(model5)

Call:
glm(formula = detectability ~ activity_cycle, data = trait_detect_clean)

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)   
(Intercept)      0.035560   0.012156   2.925  0.00991 **
activity_cycle2  0.005976   0.018378   0.325  0.74925   
activity_cycle3 -0.005584   0.024312  -0.230  0.82123   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.00132988)

    Null deviance: 0.021589  on 18  degrees of freedom
Residual deviance: 0.021278  on 16  degrees of freedom
AIC: -67.176

Number of Fisher Scoring iterations: 2
aov1 <- aov(detectability ~ activity_cycle, data = trait_detect_clean)
summary(aov1)
               Df  Sum Sq   Mean Sq F value Pr(>F)
activity_cycle  2 0.00031 0.0001552   0.117  0.891
Residuals      16 0.02128 0.0013299               
model6 <- glm(detectability ~ terrestrial_volant, data = trait_detect_clean)
summary(model6)

Call:
glm(formula = detectability ~ terrestrial_volant, data = trait_detect_clean)

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         0.036544   0.009199   3.972 0.000984 ***
terrestrial_volant1 0.001597   0.020050   0.080 0.937460    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.001269441)

    Null deviance: 0.021589  on 18  degrees of freedom
Residual deviance: 0.021580  on 17  degrees of freedom
AIC: -68.908

Number of Fisher Scoring iterations: 2

Species Richness

Given the inherent inability of passive acoustic monitoring to detect non-vocalising species, the mammal community was considered in two ways: i) entire mammal community (All Mammals), and ii) vocal mammals only (Vocal Mammals).

All Mammals

Considering the entire mammal community with vocal and non-vocal species.

Bayesian Analysis

Priors

Defining priors for Bayesian models.

Fit model

Compare Models

Diagnostics

Investigation

Methods pairwise comparison across sites
(diff.methods.avg <- methods.brm1b %>%
  emmeans(~assessment.method) %>%
  regrid(transform = "none") %>%
  pairs() %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n())) #to see if there is any change
# A tibble: 6 × 5
  contrast             `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
  <chr>                                      <dbl>       <dbl>       <dbl>                       <dbl>
1 OBM - PAM.survey                           2.46        1.85        3.12                       1     
2 OBM - camera                               1.44        1.16        1.77                       0.999 
3 PAM.all - OBM                              0.863       0.691       1.04                       0.0729
4 PAM.all - PAM.survey                       2.13        1.62        2.72                       1     
5 PAM.all - camera                           1.25        0.992       1.52                       0.977 
6 PAM.survey - camera                        0.587       0.430       0.744                      0     
(diff.methods.avg.rev <- methods.brm1b %>%
  emmeans(~assessment.method) %>%
  regrid(transform = "none") %>%
  pairs(reverse = TRUE) %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n())) #to see if there is any change
# A tibble: 6 × 5
  contrast             `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
  <chr>                                      <dbl>       <dbl>       <dbl>                       <dbl>
1 OBM - PAM.all                              1.16        0.932       1.39                     0.927   
2 PAM.survey - OBM                           0.406       0.304       0.513                    0       
3 PAM.survey - PAM.all                       0.469       0.359       0.601                    0       
4 camera - OBM                               0.693       0.555       0.848                    0.000833
5 camera - PAM.all                           0.801       0.651       0.999                    0.0229  
6 camera - PAM.survey                        1.70        1.26        2.17                     1       

Visualisation

pal <- c("#189F9F","#FF8C00", "#18709F", "#A034F0")
pal_dark <- after_scale(darken(pal, .1, space = "HLS"))

(methods.boxplot <- ggplot(mammal_data_summary, aes(assessment.method, richness)) +
  geom_boxplot(aes(color = assessment.method,
                   fill = after_scale(desaturate(lighten(color, .8), .4))),
               position = position_dodge2(width = 0.6),
               width = .6, outlier.shape = NA) +
  geom_point(aes(fill = assessment.method), color = "transparent",
             position = position_dodge2(width = 0.6),
             shape = 21, stroke = .4, size = 3.5, alpha = .3) +
  my.theme() +
  scale_y_continuous(limits = c(0,20),breaks = seq(0, 20, by = 2),
                     name = "Total Species Richness") +
  scale_x_discrete(name = "", labels = c("", "", "", "")) +
  scale_fill_manual(values = pal, guide = "none") +
  scale_colour_manual(values = pal_dark, name = "", labels = c("PAM (all audio)","OBM", "PAM (survey period)","Camera Trap")) +
  theme(legend.text = element_text(size = 14, margin = margin(t = 5)),
        legend.position = c(y=0.45, x=-0.15),
        legend.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(5, 10, 20, 10), units = "mm")) +
  guides(colour = guide_legend(nrow = 1)) +
    geom_signif(comparisons = list(c("PAM.all", "PAM.survey"),
                                   c("OBM", "PAM.survey"),
                                   c("OBM", "camera"),
                                   c("PAM.survey", "camera")),
                annotations = c("2.13, HDI 95%[1.62,2.72]",
                                "2.46, HDI 95%[1.85,3.12]",
                                "1.44, HDI 95%[1.16,1.77]",
                                "0.59, HDI 95%[0.43,0.74]"),
                y_position = c(16, 17.5, 19, 12)))

Vocal Mammals

Considering vocal mammals only and excluding non-vocal mammals from the analysis.

Bayesian Analysis

Priors

Defining priors for Bayesian models.

Fit model

Compare Models
Diagnostics

Investigation

Methods pairwise comparison across sites
(diff.methods.vocal.avg <- methods.vocal.brm1b %>%
  emmeans(~assessment.method) %>%
  regrid(transform = "none") %>%
  pairs() %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n())) #to see if there is any change
# A tibble: 6 × 5
  contrast             `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
  <chr>                                      <dbl>       <dbl>       <dbl>                       <dbl>
1 OBM - PAM.survey                            1.34       0.989        1.74                       0.980
2 OBM - camera                                1.46       1.05         1.87                       0.995
3 PAM.all - OBM                               1.58       1.24         1.98                       1    
4 PAM.all - PAM.survey                        2.12       1.61         2.69                       1    
5 PAM.all - camera                            2.29       1.73         2.93                       1    
6 PAM.survey - camera                         1.08       0.785        1.43                       0.697
(diff.methods.vocal.avg.rev <- methods.vocal.brm1b %>%
  emmeans(~assessment.method) %>%
  regrid(transform = "none") %>%
  pairs(reverse = TRUE) %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n())) #to see if there is any change
# A tibble: 6 × 5
  contrast             `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
  <chr>                                      <dbl>       <dbl>       <dbl>                       <dbl>
1 OBM - PAM.all                              0.633       0.498       0.799                      0     
2 PAM.survey - OBM                           0.746       0.555       0.974                      0.0204
3 PAM.survey - PAM.all                       0.472       0.357       0.600                      0     
4 camera - OBM                               0.686       0.514       0.905                      0.005 
5 camera - PAM.all                           0.436       0.325       0.550                      0     
6 camera - PAM.survey                        0.923       0.681       1.23                       0.303 

Visualisation

pal <- c("#189F9F","#FF8C00", "#18709F","#A034F0")
pal_dark <- after_scale(darken(pal, .1, space = "HLS"))

(methods.boxplot.vocal <- ggplot(mammal_data_vocal_summary, aes(assessment.method, richness)) +
  geom_boxplot(aes(color = assessment.method,
                   fill = after_scale(desaturate(lighten(color, .8), .4))),
               position = position_dodge2(width = 0.6),
               width = .6, outlier.shape = NA) +
  geom_point(aes(fill = assessment.method), color = "transparent",
             position = position_dodge2(width = 0.6),
             shape = 21, stroke = .4, size = 3.5, alpha = .3) +
  my.theme() +
  scale_y_continuous(limits = c(0,20),breaks = seq(0, 20, by = 2),
                     name = "") +
  scale_x_discrete(name = "", labels = c("", "", "", "")) +
  scale_fill_manual(values = pal, guide = "none") +
  scale_colour_manual(values = pal_dark, name = "", labels = c("PAM (all audio)","OBM", "PAM (survey period)","Camera Trap")) +
  theme(legend.text = element_text(size = 14, margin = margin(t = 5)),
        legend.position = c(y=0.45, x=-0.15),
        legend.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(5, 10, 20, 10), units = "mm")) +
  guides(colour = guide_legend(nrow = 1)) +
  geom_signif(comparisons = list(c("PAM.all", "OBM"),
                                 c("PAM.all", "PAM.survey"),
                                 c("PAM.all", "camera"),
                                 c("OBM", "camera")),
                annotations = c("1.58, HDI 95%[1.24,1.98]",
                                "2.12, HDI 95%[1.61,2.69]",
                                "2.29, HDI 95%[1.73,2.93]",
                                "1.46, HDI 95%[1.05,1.87]"),
                y_position = c(12.5, 14, 15.5, 9)))

Combined Figure
(methods.vocal.combined <- (methods.boxplot + ggtitle('All mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (methods.boxplot.vocal + ggtitle('Vocal mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   theme(legend.position = c(-0.25,-0.2), legend.text = element_text(size = 22),
         plot.margin = unit(c(5, 10, 30, 10), units = "mm")))

Correlation PAM vocal mammals & OBM mammals

Investigate a correlation between species richness of vocal mammals detected by PAM and species richness of the entire mammal community detected by OBM.

cor_mammal_data <- mammal_data_summary |> 
  filter(assessment.method %in% c("OBM", "PAM.all")) |> 
  group_by(site.plot) |> 
  pivot_wider(names_from = "assessment.method", #transformed dataset from long to wide format
              values_from = "richness",
              values_fill = list(number=0))

(corr_plot <- ggplot(cor_mammal_data, aes(PAM.all,OBM)) +
  geom_point(aes(col = site), size = 3) +
  geom_smooth(method = "lm") +
  my.theme() +
  theme(legend.position = "bottom") +
  scale_y_continuous(name = "Species Richness OBM", breaks = seq(4, 16, by = 4)) +
  scale_x_continuous(name = "Species Richness PAM", breaks = seq(2, 13, by = 2)) +
  scale_color_manual(name = "", values = c("seagreen", "steelblue","salmon2", "darkorchid1", "aquamarine3", "goldenrod1")) +
  guides(color = guide_legend(nrow = 1)) +
  annotate("text", x = 5, y = 15,
           label = "R^2 == 0.21", parse = TRUE, size = 5))
`geom_smooth()` using formula = 'y ~ x'

# Perform Pearson's correlation test
correlation_test <- cor.test(cor_mammal_data$PAM.all, cor_mammal_data$OBM, method = "pearson")

print(correlation_test)

    Pearson's product-moment correlation

data:  cor_mammal_data$PAM.all and cor_mammal_data$OBM
t = 2.4133, df = 22, p-value = 0.02458
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.06636028 0.72677443
sample estimates:
      cor 
0.4575097 
(R.squared <- 0.4575097^2)
[1] 0.2093151

Community Composition

All Mammals

Considering the entire mammal community with vocal and non-vocal species.

Prepare data

mammal_community3 <- mammal_data %>% 
  filter(!recapture %in% c("y")) %>%  # Removed recaptures
  select(site, assessment.method,common.name) %>% # Selected relevant grouping variables
  group_by(across(everything())) %>% # Grouped by all available columns
  mutate(number = n()) %>% # Created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% # Merged the duplicate rows
  pivot_wider(names_from = "common.name", # Transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% # Replaced NAs with 0
  arrange(site) %>% # Ordered by site
  mutate(assessment.method = factor(assessment.method)) # Turned asssessment methods into factors

In preparation for community analyses, only keep columns and rows with values above 0.

all <- mammal_community3 %>% select(where(~ any(. != 0))) %>% filter(rowSums(.[,3:ncol(.)]) > 0)

Transformation of the data to binary (presence/absence) for Jaccard dissimilarity.

all_jacc <- all %>%
  select(where(is.numeric)) %>% # Removed columns with non-numerical values
  decostand(method="pa", MARGIN = 1) %>% # Standardised presence/absence (method = "pa") by rows (MARGIN = 1)
  cbind(all[,1:2]) %>% 
  select(site, assessment.method, everything())

NMDS Ordination

Jaccard dissimilarity

set.seed(11)
all_nmds_jacc <-  metaMDS(all_jacc[4:ncol(all_jacc)], distance="jaccard", k=3,trymax=100, autotransform = FALSE)
Run 0 stress 0.1055242 
Run 1 stress 0.1207629 
Run 2 stress 0.1047353 
... New best solution
... Procrustes: rmse 0.02496561  max resid 0.0873796 
Run 3 stress 0.105525 
Run 4 stress 0.1083786 
Run 5 stress 0.109782 
Run 6 stress 0.104735 
... New best solution
... Procrustes: rmse 0.0002099691  max resid 0.0007380798 
... Similar to previous best
Run 7 stress 0.1047354 
... Procrustes: rmse 0.0002734112  max resid 0.0009616091 
... Similar to previous best
Run 8 stress 0.114386 
Run 9 stress 0.1207632 
Run 10 stress 0.1133275 
Run 11 stress 0.1097819 
Run 12 stress 0.1123875 
Run 13 stress 0.1133275 
Run 14 stress 0.1145607 
Run 15 stress 0.1133275 
Run 16 stress 0.1055252 
Run 17 stress 0.1166412 
Run 18 stress 0.1145605 
Run 19 stress 0.1155145 
Run 20 stress 0.1200923 
*** Best solution repeated 2 times

Visualisation

Jaccard dissimilarity

Function to create nMDS plots with centroids.

create_plot <- function(data, nmds_data) {
  methods <- nmds_data$points
  species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
  df <- cbind.data.frame(data[,1:2], methods)
  gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
    mutate(assessment.method = factor(assessment.method, levels = c("PAM.all", "OBM","PAM.survey", "camera")))
  ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
    geom_point(size=3) + # Smaller outside points
    geom_point(aes(x=mean.x,y=mean.y),size=5) + # This controls the centroids
    scale_colour_manual(values = c("#189F9F","#FF8C00", "#18709F","#A034F0"),
                        name = "",
                        labels = c("PAM (all audio)", "OBM", "PAM (survey period)","Camera Trap")) +
    geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) + # This controls the connections between points and centroids
    my.theme() +
    theme(legend.position = "bottom")
}

(nmds.all.plot.jacc.centroid <- create_plot(all_jacc, all_nmds_jacc))

Investigation

Adonis

Jaccard dissimilarity
all.dist2 <- vegdist(all_jacc[3:ncol(all_jacc)], 'jaccard')
(all.adonis <- adonis2(all.dist2 ~ assessment.method, data=all_jacc))
Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 999

adonis2(formula = all.dist2 ~ assessment.method, data = all_jacc)
         Df SumOfSqs      R2     F Pr(>F)    
Model     3   1.9870 0.30827 2.971  0.001 ***
Residual 20   4.4586 0.69173                 
Total    23   6.4456 1.00000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Permutation test: significant difference between assessment methods for mammal communities.

#Let's which method is different from which
(adonis.methods.jacc <- adonis.pair(all.dist2, all_jacc$assessment.method, corr.method = "holm"))
             combination SumsOfSqs   MeanSqs  F.Model        R2     P.value P.value.corrected
1         camera <-> OBM 0.4559165 0.4559165 1.692142 0.1447247 0.105894106       0.211788212
2     camera <-> PAM.all 0.9499132 0.9499132 4.958432 0.3314807 0.000999001       0.005994006
3  camera <-> PAM.survey 1.1114297 1.1114297 4.864324 0.3272483 0.003996004       0.015984016
4        OBM <-> PAM.all 0.5417048 0.5417048 2.492023 0.1994891 0.018981019       0.056943057
5     OBM <-> PAM.survey 0.6990302 0.6990302 2.748990 0.2156241 0.002997003       0.014985015
6 PAM.all <-> PAM.survey 0.2159058 0.2159058 1.223748 0.1090320 0.290709291       0.290709291
#comparing all methods with each other

Vocal Mammals

Considering vocal mammals only and excluding non-vocal mammals from the analysis.

Prepare data

mammal_community4 <- mammal_data_vocal %>% 
  filter(!recapture %in% c("y")) %>%  
  select(site, assessment.method,common.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site) %>% 
  mutate(assessment.method = factor(assessment.method))
all_vocal <- mammal_community4 %>% select(where(~ any(. != 0))) %>% filter(rowSums(.[,3:ncol(.)]) > 0)
#Function decostand() standardises presence/absence (method = "pa") by rows (MARGIN = 1)
all_vocal_jacc <- all_vocal %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1) %>% 
  cbind(all_vocal[,1:2]) %>% 
  select(site, assessment.method, everything())

NMDS Ordination

Jaccard dissimilarity

set.seed(11)
all_nmds_jacc_vocal <-  metaMDS(all_vocal_jacc[4:ncol(all_vocal_jacc)], distance="jaccard", k=3,trymax=100, autotransform = FALSE)

Visualisation

Jaccard dissimilarity

create_plot <- function(data, nmds_data) {
  methods <- nmds_data$points
  species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
  df <- cbind.data.frame(data[,1:2], methods)
  gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
    mutate(assessment.method = factor(assessment.method, levels = c("PAM.all", "OBM","PAM.survey", "camera")))
  ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
    geom_point(size=3) +
    geom_point(aes(x=mean.x,y=mean.y),size=5) +
    scale_colour_manual(values = c("#189F9F","#FF8C00", "#18709F","#A034F0"),
                        name = "",
                        labels = c("PAM (all audio)", "OBM", "PAM (survey period)","Camera Trap")) +
    geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) +
    my.theme() + theme(legend.position = "bottom")
}

(nmds.all.vocal.plot.jacc.centroid <- create_plot(all_vocal_jacc, all_nmds_jacc_vocal))

Investigation

Adonis

Jaccard dissimilarity
all.dist_vocal <- vegdist(all_vocal_jacc[3:ncol(all_vocal_jacc)], 'jaccard')
(all.adonis_vocal <- adonis2(all.dist_vocal ~ assessment.method, data=all_vocal_jacc))
Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 999

adonis2(formula = all.dist_vocal ~ assessment.method, data = all_vocal_jacc)
         Df SumOfSqs     R2      F Pr(>F)   
Model     3   1.3414 0.2559 2.2927   0.01 **
Residual 20   3.9006 0.7441                 
Total    23   5.2420 1.0000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Permutation test: significant difference between assessment methods for mammal communities.

#Let's which method is different from which
(adonis.methods.jacc_vocal <- adonis.pair(all.dist_vocal, all_vocal_jacc$assessment.method, corr.method = "holm"))
             combination SumsOfSqs   MeanSqs   F.Model         R2     P.value P.value.corrected
1         camera <-> OBM 0.4621827 0.4621827 2.1635228 0.17786976 0.071928072        0.28771229
2     camera <-> PAM.all 0.6947730 0.6947730 4.1757138 0.29456815 0.004995005        0.02497502
3  camera <-> PAM.survey 0.8980116 0.8980116 4.4172857 0.30638816 0.003996004        0.02397602
4        OBM <-> PAM.all 0.1200808 0.1200808 0.6429681 0.06041248 0.713286713        0.71328671
5     OBM <-> PAM.survey 0.2919053 0.2919053 1.3050666 0.11544086 0.233766234        0.70129870
6 PAM.all <-> PAM.survey 0.2159058 0.2159058 1.2237477 0.10903200 0.290709291        0.70129870
#comparing all methods with each other

Combined Figure

(nmds.jacc.combined <- (nmds.all.plot.jacc.centroid + ggtitle('All mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (nmds.all.vocal.plot.jacc.centroid + ggtitle('Vocal mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   theme(legend.position = c(-0.25,-0.2), legend.text = element_text(size = 22),
         plot.margin = unit(c(5, 10, 30, 10), units = "mm")) +
   guides(colour = guide_legend(nrow = 1)))

Survey Effort (Species accumulation curves)

All Mammals

Considering the entire mammal community with vocal and non-vocal species.

Prepare data

Merged duplicate rows for duplicate observations and transformed the dataset to a wide format.

mammal_community <- mammal_data %>% 
  select(site, SamplingDay, assessment.method,common.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site) %>% 
  filter(assessment.method != "PAM.all") %>% 
  mutate(SamplingDay = as.factor(SamplingDay))

Added a category for the total richness.

result.i <- vector("list",length(unique(mammal_community$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(mammal_community$site))){ #looped over each study site
  site.i <- mammal_community[mammal_community$site == unique(mammal_community$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
  for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
    
    
    SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] #subset to the SamplingDay
    
    meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("total"))
    
    SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
    
    colnames(total.a) <- colnames(mammal_community)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.i <- do.call("rbind.data.frame",result.i)

mammal_community2 <- rbind.data.frame(mammal_community, result.i)

Added a remote sensing category to assessment methods. Remote sensing was a combination of PAM and camera trapping.

result.j <- vector("list",length(unique(mammal_community$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(mammal_community$site))){ #looped over each study site
  site.i <- mammal_community[mammal_community$site == unique(mammal_community$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
  for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
    
    
    SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] %>% 
      filter(assessment.method %in% c("camera", "PAM.survey"))
    
    meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("remote"))
    
    SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
    
    colnames(total.a) <- colnames(mammal_community)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.j[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.j <- do.call("rbind.data.frame",result.j)

mammal_community2 <- rbind.data.frame(mammal_community2, result.j) %>% 
  drop_na()

Loop

Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.

# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 30, ncol = (ncol(mammal_community2)))) %>% 
  replace(is.na(.), 0) %>% 
  set_names(names(mammal_community2)) %>% 
  select(-(1:3))

result.i <- vector("list",length(unique(mammal_community2$site))) #created an empty list that was filled as the loop ran

set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(mammal_community2$site))){ #looped over each study site
  site.i <- mammal_community2[mammal_community2$site == unique(mammal_community2$site)[i],] #subset to the site
  
  dates.i <- ifelse(unique(site.i$site) %in% c("Tarcutta", "Undara", "Wambiana", "Rinyirru"), 28, 21) #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day 
  
  result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
  for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
    
    
    method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
    
    if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
      sampling.a <- method.a$assessment.method[1] #saved the name of the method
      
      method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
      
      if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
        method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
      
      accum.a <- specaccum(method.a,"random") #calculated data
      
      #created data frame of metadata, richness, and standard deviation
      accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
      colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
      
      accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
      accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
      
      result.a[[a]] <- accum.a}} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

Transformed resulting list to a data frame.

mammals_wide.result <- do.call("rbind.data.frame",result.i) %>% #compressed result.i list into data frame
mutate(assessment.method = factor(assessment.method, levels = c("OBM",
                                                                "camera",
                                                                "PAM.survey",
                                                                "remote",
                                                                "total")),
         site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

Visualisation

Plot of rarefaction curves for each method split by sites.

(specaccum_mammals <- ggplot(mammals_wide.result, aes(x = day, y = richness, col = assessment.method)) + 
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.3))),
              linetype = 0.1) +
  geom_vline(xintercept=c(7,14,21, 28), linetype="dotted", colour = "black") +
  geom_line(aes(group = assessment.method),linewidth = 1.3) +
  my.theme() +
  facet_wrap(~site) +
  theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
        strip.background = element_rect(fill = "lightgrey"),
        axis.title = element_text(size = 26),
        axis.text = element_text(size = 26),
        legend.text = element_text(size = 26),
        legend.position = "bottom",
        strip.text.x = element_text(size = 26)) +
  scale_y_continuous(limits = c(0,25), breaks = seq(0, 25, by = 5),
                     name = "Species Richness") +
  scale_x_continuous(limits = c(0,28), breaks = seq(0, 30, by = 7),
                     name = "Survey Effort (Days)") +
  scale_colour_manual(values = c("#FF8C00", "#A034F0", "#189F9F","#D14103","black"),
                      name = "",
                      labels = c("OBM","Camera Trap", "PAM (survey period)", "Remote Sensing",
                                 "Total Richness")))

Vocal Mammals

Considering vocal mammals only and excluding non-vocal mammals from the analysis.

Prepare data

Merged duplicate rows for duplicate observations and transformed the dataset to a wide format.

mammal_community_vocal <- mammal_data_vocal %>% 
  select(site, SamplingDay, assessment.method,common.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site) %>% 
  filter(assessment.method != "PAM.all") %>% 
  mutate(SamplingDay = as.factor(SamplingDay))

Added a total category to assessment methods.

result.i <- vector("list",length(unique(mammal_community_vocal$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(mammal_community_vocal$site))){ #looped over each study site
  site.i <- mammal_community_vocal[mammal_community_vocal$site == unique(mammal_community_vocal$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
  for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
    
    
    SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] #subset to the SamplingDay
    
    meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("total"))
    
    SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
    
    colnames(total.a) <- colnames(mammal_community_vocal)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.i <- do.call("rbind.data.frame",result.i)

mammal_community_vocal2 <- rbind.data.frame(mammal_community_vocal, result.i)

Added a remote sensing category to assessment methods. Remote sensing was a combination of PAM and camera trapping.

result.j <- vector("list",length(unique(mammal_community_vocal$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(mammal_community_vocal$site))){ #looped over each study site
  site.i <- mammal_community_vocal[mammal_community_vocal$site == unique(mammal_community_vocal$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
  for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
    
    
    SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] %>% 
      filter(assessment.method %in% c("camera", "PAM.survey"))
    
    meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("remote"))
    
    SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
    
    colnames(total.a) <- colnames(mammal_community_vocal)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.j[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.j <- do.call("rbind.data.frame",result.j)

mammal_community_vocal2 <- rbind.data.frame(mammal_community_vocal2, result.j) %>% 
  drop_na()

Loop

Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.

# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 30, ncol = (ncol(mammal_community_vocal2)))) %>% 
  replace(is.na(.), 0) %>% 
  set_names(names(mammal_community_vocal2)) %>% 
  select(-(1:3))

result.i <- vector("list",length(unique(mammal_community_vocal2$site))) #created an empty list that was filled as the loop ran

set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(mammal_community_vocal2$site))){ #looped over each study site
  site.i <- mammal_community_vocal2[mammal_community_vocal2$site == unique(mammal_community_vocal2$site)[i],] #subset to the site
  
  dates.i <- ifelse(unique(site.i$site) %in% c("Tarcutta", "Undara", "Wambiana", "Rinyirru"), 28, 21) #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day 
  
  result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
  for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
    
    
    method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
    
    if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
      sampling.a <- method.a$assessment.method[1] #saved the name of the method
      
      method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
      
      if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
        method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
      
      accum.a <- specaccum(method.a,"random") #calculated data
      
      #created data frame of metadata, richness, and standard deviation
      accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
      colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
      
      accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
      accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
      
      result.a[[a]] <- accum.a}} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

Transformed resulting list to a data frame.

mammals_wide.result_vocal <- do.call("rbind.data.frame",result.i) %>% #compressed result.i list into data frame
mutate(assessment.method = factor(assessment.method, levels = c("OBM",
                                                                "camera",
                                                                "PAM.survey",
                                                                "remote",
                                                                "total")),
         site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

Visualisation

Plot of rarefaction curves for each method split by sites.

(specaccum_mammals_vocal <- ggplot(mammals_wide.result_vocal, aes(x = day, y = richness, col = assessment.method)) + 
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.3))),
              linetype = 0.1) +
  geom_vline(xintercept=c(7,14,21, 28), linetype="dotted", colour = "black") +
  geom_line(aes(group = assessment.method), linewidth = 1.3) +
  my.theme() +
  facet_wrap(~site) +
  theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
        strip.background = element_rect(fill = "lightgrey"),
        axis.title = element_text(size = 26),
        axis.text = element_text(size = 26),
        legend.text = element_text(size = 26),
        legend.position = "bottom",
        strip.text.x = element_text(size = 26)) +
  scale_y_continuous(limits = c(0,15), breaks = seq(0, 15, by = 5),
                     name = "") +
  scale_x_continuous(limits = c(0,28), breaks = seq(0, 30, by = 7),
                     name = "Survey Effort (Days)") +
  scale_colour_manual(values = c("#FF8C00", "#A034F0", "#189F9F","#D14103" ,"black"),
                      name = "",
                      labels = c("OBM","Camera Trap", "PAM (survey period)", "Remote Sensing",
                                 "Total Richness")))

Combined Figure

(specaccum_mammals_combined <- (specaccum_mammals + ggtitle('All mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (specaccum_mammals_vocal + ggtitle('Vocal mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   theme(legend.position = c(-0.1,-0.2), legend.text = element_text(size = 22),
         plot.margin = unit(c(5, 10, 30, 10), units = "mm")) +
   guides(colour = guide_legend(nrow = 1)))

PAM survey length

How many days until PAM.all reaches values similar to OBM for all mammals and only vocal mammals. In other words how long should we leave recorders out on average until we get the best number of species.

all_mammals <- mammal_data %>% 
  filter(assessment.method == "PAM.all") %>% 
  select(site, SamplingDay, assessment.method,common.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site) %>%
  mutate(SamplingDay = as.factor(SamplingDay))

Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.

# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 1235, ncol = (ncol(all_mammals)))) %>% 
  replace(is.na(.), 0) %>% 
  set_names(names(all_mammals)) %>% 
  select(-(1:3))

result.i <- vector("list",length(unique(all_mammals$site))) #created an empty list that was filled as the loop ran

set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(all_mammals$site))){ #looped over each study site
  site.i <- all_mammals[all_mammals$site == unique(all_mammals$site)[i],] #subset to the site
  
   #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day 
  site_values <- c(Tarcutta = 986, Duval = 679, Mourachan = 783, Wambiana = 1235, Undara = 780, Rinyirru = 488)
  
  dates.i <- site_values[match(unique(site.i$site), names(site_values))]
  
  result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
  for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
    
    
    method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
    
    if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
      sampling.a <- method.a$assessment.method[1] #saved the name of the method
      
      method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
      
      if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
        method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
      
      accum.a <- specaccum(method.a,"random") #calculated data
      
      #created data frame of metadata, richness, and standard deviation
      accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
      colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
      
      accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
      accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
      
      result.a[[a]] <- accum.a}} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

Transformed resulting list to a data frame.

all_mammals.result <- do.call("rbind.data.frame",result.i) %>%  #compressed result.i list into data
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))
pam.max <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == max(day)) %>% 
  ungroup() %>% 
  select(site, richness.max = richness)

pam.28 <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == 28) %>% 
  ungroup() %>% 
  select(richness.28 = richness) %>% 
  mutate(prop.28 = round(richness.28/pam.max$richness.max, digits = 2))

pam.90 <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == 90) %>% 
  ungroup() %>% 
  select(richness.90 = richness) %>% 
  mutate(prop.90 = round(richness.90/pam.max$richness.max, digits = 2))

pam.180 <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == 180) %>% 
  ungroup() %>% 
  select(richness.180 = richness) %>% 
  mutate(prop.180 = round(richness.180/pam.max$richness.max, digits = 2))

pam.365 <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == 365) %>% 
  ungroup() %>% 
  select(richness.365 = richness) %>% 
  mutate(prop.365 = round(richness.365/pam.max$richness.max, digits = 2))

pam.rich.prop <- cbind(pam.max,pam.28,pam.90,pam.180,pam.365) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

(sampling.table <- pam.rich.prop %>% 
  select(Site = site, "PAM (total richness)" = richness.max,
         "PAM (proportion 28 days)" = prop.28,
         "PAM (proportion 90 days)" = prop.90,
         "PAM (proportion 180 days)" = prop.180,
         "PAM (proportion 365 days)" = prop.365) %>% 
  add_row(Site = "Average", "PAM (total richness)" = round(mean(pam.rich.prop$richness.max)),
          "PAM (proportion 28 days)" = round(mean(pam.rich.prop$prop.28), digits = 2),
          "PAM (proportion 90 days)" = round(mean(pam.rich.prop$prop.90), digits = 2),
          "PAM (proportion 180 days)" = round(mean(pam.rich.prop$prop.180), digits = 2),
          "PAM (proportion 365 days)" = round(mean(pam.rich.prop$prop.365), digits = 2)))
       Site PAM (total richness) PAM (proportion 28 days) PAM (proportion 90 days) PAM (proportion 180 days)
1     Duval                   10                     0.72                     0.86                      0.91
2 Mourachan                    7                     0.68                     0.86                      0.89
3  Rinyirru                    9                     0.61                     0.79                      0.87
4  Tarcutta                   13                     0.66                     0.78                      0.86
5    Undara                    7                     0.50                     0.75                      0.86
6  Wambiana                   10                     0.79                     0.91                      0.96
7   Average                    9                     0.66                     0.82                      0.89
  PAM (proportion 365 days)
1                      0.95
2                      0.92
3                      0.96
4                      0.94
5                      0.92
6                      0.99
7                      0.95

Visualisation

vline_data <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
                         xintercept = c(12, 175, 17, 88, 676,320)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

vline_data2 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
                         xintercept = c(12, 35, 12, 17, 12,9)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

hline_data_vocal <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(4.1, 6, 7, 6,12,10)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

hline_data_vocal2 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(3.9, 4, 6, 4,6,6)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))



hline_data_cam <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(4, 4, 6, 4,6,6)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

vline_data_cam <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
                         xintercept = c(12, 35, 12, 17, 12,9)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

hline_data_28 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(5.48, 3.51, 7.88, 4.32,8.52,6.87)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

vline_data_28 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
                         xintercept = c(28, 28, 28, 21, 21,28)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

hline_data_total <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(8, 7, 11, 8,13,11)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

vline_data_total <- data.frame(site = c("Rinyirru", "Undara","Tarcutta"),
                         xintercept = c(210, 771, 968)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara","Tarcutta")))


(spec_all_mammals_plot <- ggplot(all_mammals.result, aes(x = day, y = richness, col = assessment.method)) +
    geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.5))),
              linetype = 0.1) +
    geom_line(aes(group = assessment.method), linewidth = 1.5) +
    my.theme() +
    facet_wrap(~site, scales = "free_x") +
    geom_hline(data = hline_data_vocal, aes(yintercept = yintercept), linetype="dashed", color = "#FF8C00",
               linewidth = 1) +
    geom_hline(data = hline_data_vocal2, aes(yintercept = yintercept), linetype="dashed", color = "#A034F0",
               linewidth = 1) +
    #geom_vline(data = vline_data, aes(xintercept = xintercept), linetype="dotted", color = "#FF8C00") +
    #geom_vline(data = vline_data2, aes(xintercept = xintercept), linetype="dotted", color = "#A034F0") +
    geom_text(data = vline_data, aes(x = 100, y = 14, label = paste0("DTM: ",xintercept)),
                                     vjust = -1, size = 8,inherit.aes = FALSE, color = "#FF8C00") +
    geom_text(data = vline_data2, aes(x = 100, y = 13, label = paste0("DTM: ",xintercept)),
                                     vjust = -1, size = 8,inherit.aes = FALSE, color = "#A034F0") +
    theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
        strip.background = element_rect(fill = "lightgrey"),
        axis.title = element_text(size = 26),
        axis.text = element_text(size = 26),
        legend.text = element_text(size = 26),
        legend.position = "bottom",
        strip.text.x = element_text(size = 26)) + 
    scale_y_continuous(limits = c(0,15), breaks = seq(0, 15, by = 5),
                     name = "Species Richness") +
    scale_x_continuous(name = "Survey Effort (Days)") +
    scale_colour_manual(values = c("#189F9F"),
                      name = "",
                      labels = c("PAM (all audio)")))

Session Info

sessionInfo()
R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Brisbane
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] pilot_4.0.1         paletteer_1.6.0     cowplot_1.1.3       ggimage_0.3.3       ggsignif_0.6.4     
 [6] report_0.6.1        HDInterval_0.2.4    patchwork_1.3.0     EcolUtils_0.1       vegan_2.6-10       
[11] lattice_0.22-6      permute_0.9-7       tidybayes_3.0.7     emmeans_1.11.1      DHARMa_0.4.7       
[16] ggeffects_2.2.1     rstan_2.32.7        StanHeaders_2.32.10 brms_2.22.0         Rcpp_1.0.14        
[21] cmdstanr_0.9.0      colorspace_2.1-1    lubridate_1.9.4     forcats_1.0.0       stringr_1.5.1      
[26] dplyr_1.1.4         purrr_1.0.4         readr_2.1.5         tidyr_1.3.1         tibble_3.2.1       
[31] ggplot2_3.5.2       tidyverse_2.0.0    

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3   tensorA_0.36.2.1     rstudioapi_0.17.1    jsonlite_2.0.0       magrittr_2.0.3      
  [6] estimability_1.5.1   magick_2.8.6         farver_2.1.2         nloptr_2.2.1         rmarkdown_2.29      
 [11] fs_1.6.6             vctrs_0.6.5          minqa_1.2.8          htmltools_0.5.8.1    distributional_0.5.0
 [16] curl_6.2.2           gridGraphics_0.5-1   htmlwidgets_1.6.4    plyr_1.8.9           mime_0.13           
 [21] iterators_1.0.14     lifecycle_1.0.4      pkgconfig_2.0.3      gap_1.6              Matrix_1.7-3        
 [26] R6_2.6.1             fastmap_1.2.0        shiny_1.10.0         rbibutils_2.3        digest_0.6.37       
 [31] rematch2_2.1.2       ps_1.9.1             qgam_2.0.0           labeling_0.4.3       timechange_0.3.0    
 [36] abind_1.4-8          mgcv_1.9-1           compiler_4.5.0       doParallel_1.0.17    bit64_4.6.0-1       
 [41] withr_3.0.2          backports_1.5.0      inline_0.3.21        QuickJSR_1.7.0       pkgbuild_1.4.7      
 [46] MASS_7.3-65          loo_2.8.0            tools_4.5.0          httpuv_1.6.16        glue_1.8.0          
 [51] promises_1.3.2       nlme_3.1-168         grid_4.5.0           checkmate_2.3.2      cluster_2.1.8.1     
 [56] reshape2_1.4.4       generics_0.1.4       gtable_0.3.6         tzdb_0.5.0           data.table_1.17.2   
 [61] hms_1.1.3            utf8_1.2.5           foreach_1.5.2        pillar_1.10.2        ggdist_3.3.3        
 [66] yulab.utils_0.2.0    vroom_1.6.5          later_1.4.2          posterior_1.6.1      splines_4.5.0       
 [71] bit_4.6.0            tidyselect_1.2.1     knitr_1.50           reformulas_0.4.1     arrayhelpers_1.1-0  
 [76] gridExtra_2.3        V8_6.0.5             stats4_4.5.0         xfun_0.52            bridgesampling_1.1-2
 [81] matrixStats_1.5.0    stringi_1.8.7        ggfun_0.1.8          yaml_2.3.10          boot_1.3-31         
 [86] evaluate_1.0.3       codetools_0.2-20     ggplotify_0.1.2      cli_3.6.5            RcppParallel_5.1.10 
 [91] xtable_1.8-4         Rdpack_2.6.4         processx_3.8.6       coda_0.19-4.1        svUnit_1.0.6        
 [96] parallel_4.5.0       rstantools_2.4.0     gap.datasets_0.0.6   bayesplot_1.12.0     Brobdingnag_1.2-9   
[101] lme4_1.1-37          mvtnorm_1.3-3        scales_1.4.0         insight_1.3.1        crayon_1.5.3        
[106] rlang_1.1.6         

Cite Packages

cite_packages()
  - Bürkner P (2017). "brms: An R Package for Bayesian Multilevel Models Using Stan." _Journal of Statistical Software_, *80*(1), 1-28. doi:10.18637/jss.v080.i01 <https://doi.org/10.18637/jss.v080.i01>. Bürkner P (2018). "Advanced Bayesian Multilevel Modeling with the R Package brms." _The R Journal_, *10*(1), 395-411. doi:10.32614/RJ-2018-017 <https://doi.org/10.32614/RJ-2018-017>. Bürkner P (2021). "Bayesian Item Response Modeling in R with brms and Stan." _Journal of Statistical Software_, *100*(5), 1-54. doi:10.18637/jss.v100.i05 <https://doi.org/10.18637/jss.v100.i05>.
  - Constantin A, Patil I (2021). "ggsignif: R Package for Displaying Significance Brackets for 'ggplot2'." _PsyArxiv_. doi:10.31234/osf.io/7awm6 <https://doi.org/10.31234/osf.io/7awm6>, <https://psyarxiv.com/7awm6>.
  - Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2025). _Rcpp: Seamless R and C++ Integration_. doi:10.32614/CRAN.package.Rcpp <https://doi.org/10.32614/CRAN.package.Rcpp>, R package version 1.0.14, <https://CRAN.R-project.org/package=Rcpp>. Eddelbuettel D, François R (2011). "Rcpp: Seamless R and C++ Integration." _Journal of Statistical Software_, *40*(8), 1-18. doi:10.18637/jss.v040.i08 <https://doi.org/10.18637/jss.v040.i08>. Eddelbuettel D (2013). _Seamless R and C++ Integration with Rcpp_. Springer, New York. doi:10.1007/978-1-4614-6868-4 <https://doi.org/10.1007/978-1-4614-6868-4>, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). "Extending R with C++: A Brief Introduction to Rcpp." _The American Statistician_, *72*(1), 28-36. doi:10.1080/00031305.2017.1375990 <https://doi.org/10.1080/00031305.2017.1375990>.
  - Gabry J, Češnovar R, Johnson A, Bronder S (2025). _cmdstanr: R Interface to 'CmdStan'_. R package version 0.9.0, <https://mc-stan.org/cmdstanr/>.
  - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate." _Journal of Statistical Software_, *40*(3), 1-25. <https://www.jstatsoft.org/v40/i03/>.
  - Hartig F (2024). _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models_. doi:10.32614/CRAN.package.DHARMa <https://doi.org/10.32614/CRAN.package.DHARMa>, R package version 0.4.7, <https://CRAN.R-project.org/package=DHARMa>.
  - Hawkins O (2025). _pilot: A minimal ggplot2 theme with an accessible discrete color palette_. R package version 4.0.1, commit 532e68c5135fc847130a76547f3d32832bf43f08, <https://github.com/olihawkins/pilot>.
  - Hvitfeldt E (2021). _paletteer: Comprehensive Collection of Color Palettes_. R package version 1.3.0, <https://github.com/EmilHvitfeldt/paletteer>.
  - Kay M (2024). _tidybayes: Tidy Data and Geoms for Bayesian Models_. doi:10.5281/zenodo.1308151 <https://doi.org/10.5281/zenodo.1308151>, R package version 3.0.7, <http://mjskay.github.io/tidybayes/>.
  - Lenth R (2025). _emmeans: Estimated Marginal Means, aka Least-Squares Means_. doi:10.32614/CRAN.package.emmeans <https://doi.org/10.32614/CRAN.package.emmeans>, R package version 1.11.1, <https://CRAN.R-project.org/package=emmeans>.
  - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects from Regression Models." _Journal of Open Source Software_, *3*(26), 772. doi:10.21105/joss.00772 <https://doi.org/10.21105/joss.00772>.
  - Makowski D, Lüdecke D, Patil I, Thériault R, Ben-Shachar M, Wiernik B (2023). "Automated Results Reporting as a Practical Tool to Improve Reproducibility and Methodological Best Practices Adoption." _CRAN_. <https://easystats.github.io/report/>.
  - Meredith M, Kruschke J (2022). _HDInterval: Highest (Posterior) Density Intervals_. doi:10.32614/CRAN.package.HDInterval <https://doi.org/10.32614/CRAN.package.HDInterval>, R package version 0.2.4, <https://CRAN.R-project.org/package=HDInterval>.
  - Müller K, Wickham H (2023). _tibble: Simple Data Frames_. doi:10.32614/CRAN.package.tibble <https://doi.org/10.32614/CRAN.package.tibble>, R package version 3.2.1, <https://CRAN.R-project.org/package=tibble>.
  - Oksanen J, Simpson G, Blanchet F, Kindt R, Legendre P, Minchin P, O'Hara R, Solymos P, Stevens M, Szoecs E, Wagner H, Barbour M, Bedward M, Bolker B, Borcard D, Carvalho G, Chirico M, De Caceres M, Durand S, Evangelista H, FitzJohn R, Friendly M, Furneaux B, Hannigan G, Hill M, Lahti L, McGlinn D, Ouellette M, Ribeiro Cunha E, Smith T, Stier A, Ter Braak C, Weedon J, Borman T (2025). _vegan: Community Ecology Package_. doi:10.32614/CRAN.package.vegan <https://doi.org/10.32614/CRAN.package.vegan>, R package version 2.6-10, <https://CRAN.R-project.org/package=vegan>.
  - Pedersen T (2024). _patchwork: The Composer of Plots_. doi:10.32614/CRAN.package.patchwork <https://doi.org/10.32614/CRAN.package.patchwork>, R package version 1.3.0, <https://CRAN.R-project.org/package=patchwork>.
  - R Core Team (2025). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.
  - Salazar G (2025). _EcolUtils: Utilities for community ecology analysis_. R package version 0.1, commit 2b57a8a8b76eba0852c56ce30f93d92c30320161, <https://github.com/GuillemSalazar/EcolUtils>.
  - Sarkar D (2008). _Lattice: Multivariate Data Visualization with R_. Springer, New York. ISBN 978-0-387-75968-5, <http://lmdvr.r-forge.r-project.org>.
  - Simpson G (2022). _permute: Functions for Generating Restricted Permutations of Data_. doi:10.32614/CRAN.package.permute <https://doi.org/10.32614/CRAN.package.permute>, R package version 0.9-7, <https://CRAN.R-project.org/package=permute>.
  - Stan Development Team (2020). "StanHeaders: Headers for the R interface to Stan." R package version 2.21.0-6, <https://mc-stan.org/>.
  - Stan Development Team (2025). "RStan: the R interface to Stan." R package version 2.32.7, <https://mc-stan.org/>.
  - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_. Springer-Verlag New York. ISBN 978-3-319-24277-4, <https://ggplot2.tidyverse.org>.
  - Wickham H (2023). _forcats: Tools for Working with Categorical Variables (Factors)_. doi:10.32614/CRAN.package.forcats <https://doi.org/10.32614/CRAN.package.forcats>, R package version 1.0.0, <https://CRAN.R-project.org/package=forcats>.
  - Wickham H (2023). _stringr: Simple, Consistent Wrappers for Common String Operations_. doi:10.32614/CRAN.package.stringr <https://doi.org/10.32614/CRAN.package.stringr>, R package version 1.5.1, <https://CRAN.R-project.org/package=stringr>.
  - Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." _Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
  - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar of Data Manipulation_. doi:10.32614/CRAN.package.dplyr <https://doi.org/10.32614/CRAN.package.dplyr>, R package version 1.1.4, <https://CRAN.R-project.org/package=dplyr>.
  - Wickham H, Henry L (2025). _purrr: Functional Programming Tools_. doi:10.32614/CRAN.package.purrr <https://doi.org/10.32614/CRAN.package.purrr>, R package version 1.0.4, <https://CRAN.R-project.org/package=purrr>.
  - Wickham H, Hester J, Bryan J (2024). _readr: Read Rectangular Text Data_. doi:10.32614/CRAN.package.readr <https://doi.org/10.32614/CRAN.package.readr>, R package version 2.1.5, <https://CRAN.R-project.org/package=readr>.
  - Wickham H, Vaughan D, Girlich M (2024). _tidyr: Tidy Messy Data_. doi:10.32614/CRAN.package.tidyr <https://doi.org/10.32614/CRAN.package.tidyr>, R package version 1.3.1, <https://CRAN.R-project.org/package=tidyr>.
  - Wilke C (2024). _cowplot: Streamlined Plot Theme and Plot Annotations for 'ggplot2'_. doi:10.32614/CRAN.package.cowplot <https://doi.org/10.32614/CRAN.package.cowplot>, R package version 1.1.3, <https://CRAN.R-project.org/package=cowplot>.
  - Yu G (2023). _ggimage: Use Image in 'ggplot2'_. doi:10.32614/CRAN.package.ggimage <https://doi.org/10.32614/CRAN.package.ggimage>, R package version 0.3.3, <https://CRAN.R-project.org/package=ggimage>.
  - Zeileis A, Fisher JC, Hornik K, Ihaka R, McWhite CD, Murrell P, Stauffer R, Wilke CO (2020). "colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes." _Journal of Statistical Software_, *96*(1), 1-49. doi:10.18637/jss.v096.i01 <https://doi.org/10.18637/jss.v096.i01>. Zeileis A, Hornik K, Murrell P (2009). "Escaping RGBland: Selecting Colors for Statistical Graphics." _Computational Statistics & Data Analysis_, *53*(9), 3259-3270. doi:10.1016/j.csda.2008.11.033 <https://doi.org/10.1016/j.csda.2008.11.033>. Stauffer R, Mayr GJ, Dabernig M, Zeileis A (2009). "Somewhere over the Rainbow: How to Make Effective Use of Colors in Meteorological Visualizations." _Bulletin of the American Meteorological Society_, *96*(2), 203-216. doi:10.1175/BAMS-D-13-00155.1 <https://doi.org/10.1175/BAMS-D-13-00155.1>.